Intro on the history of Charlottesville… the housing market… renter explotation…
Loading packages for exploration
library(readr)
library(leaflet)
library(ggplot2)
library(GGally)
## Warning: package 'GGally' was built under R version 4.2.1
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(scales)
##
## Attaching package: 'scales'
## The following object is masked from 'package:readr':
##
## col_factor
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✔ tibble 3.1.7 ✔ stringr 1.4.0
## ✔ tidyr 1.2.0 ✔ forcats 0.5.1
## ✔ purrr 0.3.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ scales::col_factor() masks readr::col_factor()
## ✖ purrr::discard() masks scales::discard()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(tigris)
## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1; sf_use_s2() is TRUE
library(viridis)
## Loading required package: viridisLite
##
## Attaching package: 'viridis'
## The following object is masked from 'package:scales':
##
## viridis_pal
library(RColorBrewer)
# Reading in exported data to R markdown
library(readr)
ACSData2020 <- read_csv("ACSData2020.csv") %>%
mutate(GEOID=as.character(GEOID))
## Rows: 317 Columns: 29
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): NAME, FipCodes, County, Region
## dbl (23): GEOID, MedRentE, MedRentM, MedTaxE, MedTaxM, RentTaxRatE, RentTaxR...
## lgl (2): PacIslMedHHIncomeE, PacIslMedHHIncomeM
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
ACSData2020
## # A tibble: 317 × 29
## GEOID NAME FipCodes County Region MedRentE MedRentM MedTaxE MedTaxM
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 51540000201 Census … 540 Charl… Charl… 1444 111 1946 655
## 2 51540000202 Census … 540 Charl… Charl… 1228 187 1826 129
## 3 51540000302 Census … 540 Charl… Charl… 1201 140 3228 368
## 4 51540000401 Census … 540 Charl… Charl… 679 115 2526 305
## 5 51540000402 Census … 540 Charl… Charl… 965 110 2364 398
## 6 51540000501 Census … 540 Charl… Charl… 1160 100 1858 219
## 7 51540000502 Census … 540 Charl… Charl… 1285 121 2812 402
## 8 51540000600 Census … 540 Charl… Charl… 1161 54 1985 619
## 9 51540000700 Census … 540 Charl… Charl… 994 133 4673 532
## 10 51540000800 Census … 540 Charl… Charl… 911 151 3184 211
## # … with 307 more rows, and 20 more variables: RentTaxRatE <dbl>,
## # RentTaxRatMOE <dbl>, MedRenterHHIncomeE <dbl>, MedRenterHHIncomeM <dbl>,
## # MedRentBurdenPercE <dbl>, MedRentBurdenPercM <dbl>, MedHHIncomeE <dbl>,
## # MedHHIncomeM <dbl>, WhiteMedHHIncomeE <dbl>, WhiteMedHHIncomeM <dbl>,
## # BlackMedHHIncomeE <dbl>, BlackMedHHIncomeM <dbl>,
## # NatAmerMedHHIncomeE <dbl>, NatAmerMedHHIncomeM <dbl>,
## # AsianMedHHIncomeE <dbl>, AsianMedHHIncomeM <dbl>, …
Let us find the median values per region.
ACSData2020 %>%
group_by(Region) %>%
summarize(MedRentEst= median(MedRentE, na.rm=TRUE), MedTaxEst=median(MedTaxE, na.rm=TRUE), MedRentTaxRatEst=median(RentTaxRatE, na.rm=TRUE), MedRenterHHIncomeEst =median(MedRenterHHIncomeE, na.rm=TRUE), MedHHIncomeEst=median(MedHHIncomeE, na.rm=TRUE), MedRentBurdenPercEst=median(MedRentBurdenPercE, na.rm=TRUE), WhiteMedHHIncomeEst=median(WhiteMedHHIncomeE, na.rm=TRUE), BlackMedHHIncomeEst=median(BlackMedHHIncomeE, na.rm=TRUE), HispanicLatinoMedHHIncomeEst=median(HispanicLatinoMedHHIncomeE, na.rm=TRUE))
## # A tibble: 2 × 10
## Region MedRentEst MedTaxEst MedRentTaxRatEst MedRenterHHInco… MedHHIncomeEst
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Charlot… 1223 2218 0.559 49860 74480
## 2 Richmond 1174 1875 0.660 48382 65958
## # … with 4 more variables: MedRentBurdenPercEst <dbl>,
## # WhiteMedHHIncomeEst <dbl>, BlackMedHHIncomeEst <dbl>,
## # HispanicLatinoMedHHIncomeEst <dbl>
#med rent, med rent tax, med rent burden, gini index, total above pov, total below pov
visulizing some of these findings…
ggplot(ACSData2020) +
aes(x = Region, y = RentTaxRatE) +
geom_boxplot()
## Warning: Removed 23 rows containing non-finite values (stat_boxplot).
ggplot(ACSData2020) +
aes(x = Region, y = MedRentBurdenPercE) +
geom_boxplot()
## Warning: Removed 13 rows containing non-finite values (stat_boxplot).
##{}
ggplot(ACSData2020, aes(RentTaxRatE))+
geom_histogram()+
facet_wrap(~Region)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 23 rows containing non-finite values (stat_bin).
ggplot(ACSData2020, aes(MedRentBurdenPercE))+
geom_histogram()+
facet_wrap(~Region)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 13 rows containing non-finite values (stat_bin).
ggplot(ACSData2020)+
aes(x=MedRentBurdenPercE, y=RentTaxRatE ) +
geom_point()+
geom_smooth()+
facet_wrap(~Region)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 23 rows containing non-finite values (stat_smooth).
## Warning: Removed 23 rows containing missing values (geom_point).
Comparing rent exploitation and rent burden to see if there is a common
problem such as those that are greater exploited also being greater
burdened.
ggplot(ACSData2020) +
aes(x = MedRentBurdenPercE) +
geom_bar()
## Warning: Removed 13 rows containing non-finite values (stat_count).
… and further by county.
ACSData2020 %>%
group_by(County) %>%
summarize(MedRentEst= median(MedRentE, na.rm=TRUE), MedTaxEst=median(MedTaxE, na.rm=TRUE), MedRentTaxRatEst=median(RentTaxRatE, na.rm=TRUE), MedRenterHHIncomeEst =median(MedRenterHHIncomeE, na.rm=TRUE), MedHHIncomeEst=median(MedHHIncomeE, na.rm=TRUE), MedRentBurdenPercEst=median(MedRentBurdenPercE, na.rm=TRUE), WhiteMedHHIncomeEst=median(WhiteMedHHIncomeE, na.rm=TRUE), BlackMedHHIncomeEst=median(BlackMedHHIncomeE, na.rm=TRUE), HispanicLatinoMedHHIncomeEst=median(HispanicLatinoMedHHIncomeE, na.rm=TRUE))
## # A tibble: 11 × 10
## County MedRentEst MedTaxEst MedRentTaxRatEst MedRenterHHInco… MedHHIncomeEst
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Albema… 1323 2657 0.542 55964 87000
## 2 Charlo… 1181 2622. 0.433 40104 62478.
## 3 Cheste… 1333 1935 0.719 60040 80102
## 4 Fluvan… 1419 1786 0.694 49581 79590
## 5 Greene 974. 1636. 0.620 51240. 64824.
## 6 Henric… 1214 1923 0.659 53133 69827
## 7 Hopewe… 910. 1104. 0.845 28625 42282
## 8 Louisa 875 1432 0.644 46964 61842
## 9 Nelson 918. 1357 0.710 44754 53579
## 10 Peters… 952 1067 0.769 34167 40552
## 11 Richmo… 1085 2148 0.484 37975 53216.
## # … with 4 more variables: MedRentBurdenPercEst <dbl>,
## # WhiteMedHHIncomeEst <dbl>, BlackMedHHIncomeEst <dbl>,
## # HispanicLatinoMedHHIncomeEst <dbl>
visualizing some of these findings…
ggplot(ACSData2020) +
aes(x = County, y = RentTaxRatE) +
geom_boxplot()
## Warning: Removed 23 rows containing non-finite values (stat_boxplot).
ggplot(ACSData2020) +
aes(x = County, y = MedRentBurdenPercE) +
geom_boxplot()
## Warning: Removed 13 rows containing non-finite values (stat_boxplot).
##{}
ggplot(ACSData2020, aes(RentTaxRatE))+
geom_histogram()+
facet_wrap(~County)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 23 rows containing non-finite values (stat_bin).
ggplot(ACSData2020, aes(MedRentBurdenPercE))+
geom_histogram()+
facet_wrap(~County)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 13 rows containing non-finite values (stat_bin).
ggplot(ACSData2020)+
aes(x=MedRentBurdenPercE, y=RentTaxRatE ) +
geom_point()+
geom_smooth()+
facet_wrap(~County)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 23 rows containing non-finite values (stat_smooth).
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 18.137
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 3.1635
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 891.83
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small. fewer
## data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 18.137
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 3.1635
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 891.83
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 17.747
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 9.553
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 3.073
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small. fewer
## data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 17.747
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 9.553
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 3.073
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in sqrt(sum.squares/one.delta): NaNs produced
## Warning in stats::qt(level/2 + 0.5, pred$df): NaNs produced
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 13.513
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 14.287
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 978.85
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small. fewer
## data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 13.513
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 14.287
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 978.85
## Warning: Removed 23 rows containing missing values (geom_point).
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
counties<-c("Albemarle", "Charlottesville", "Fluvanna", "Greene", "Louisa","Nelson", "Chesterfield County", "Henrico County", "Hopewell City", "Petersburg City", "Richmond City")
countytracts <- tracts(state = "51", county = counties, year = 2020)
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%
As defined by Crowell, rent exploitation is the ratio of median rents to median property taxes (Crowell 2). We will examine how rent exploitation varies by counties in the Charlottesville and Richmond regions.
ACSData2020<-countytracts %>%
left_join(ACSData2020, by="GEOID")
ACSData2020 <- st_transform(ACSData2020, 4326)
pal <- colorNumeric("viridis", reverse = FALSE, domain = ACSData2020$RentTaxRatE) # viridis
leaflet(ACSData2020) %>%
#addTiles() %>% # basemap
addProviderTiles("CartoDB.Positron") %>% # CartoDB.Positron, Stamen.Toner, Esri.NatGeoWorldMap
addPolygons(fillColor = ~pal(RentTaxRatE),
fillOpacity = 0.8,
color = "white",
weight = 2,
highlight = highlightOptions(
weight = 3,
fillOpacity = 0.9,
bringToFront = T),
popup = paste0("Region: ", ACSData2020$Region, "<br>",
"County: ", ACSData2020$County, "<br>",
"Tract: ", ACSData2020$NAMELSAD, "<br>",
"Ratio: ", round(ACSData2020$RentTaxRatE, 2)))%>%
addLegend("topright", pal = pal, values = ACSData2020$RentTaxRatE,
title = "Rent Exploitation", opacity = 0.7)
Recalling that if rent exploitation> 1 then the tenant is paying more than the property is worth judging by how much the owner is paying for its property tax.
#create a filtered map that explicity shows those financially exploited?
ACSData2020 %>%
ggplot()+
geom_sf(aes(fill= RentTaxRatE))+
#geom_sf_label(aes(label = County), size=3, color ="black", fontface = "bold")+
theme_void()
ACSData2020 %>%
#filter(Region == "Charlottesville") %>%
ggplot()+
geom_sf(aes(fill=RentTaxRatE))+
scale_fill_viridis_c()+
scale_color_viridis_c()+
ggtitle("Rent Exploitation in Charlottesville and Richmond, VA", subtitle= "Data from the 2020 American Community Survey (ACS)")+
labs(fill="Rent Exploitation")
ggsave("Rent Exploitation Charlottesville and Richmond.png")
## Saving 7 x 5 in image
ACSData2020 %>%
filter(Region == "Richmond") %>%
ggplot()+
geom_sf(aes(fill=RentTaxRatE))+
scale_fill_viridis_c()+
scale_color_viridis_c()+
ggtitle("Rent Explotation in Richmond, VA", subtitle= "Data from the 2020 American Community Survey")+
labs(fill="Rent Exploitation")
Additionally, rent burden is the ratio of median renters income to
median rent as a measurement of rent affordabilty as if the renter is
paying greater than 30% of their income they range from burdened to
extremely burdened in terms of being able to spend their income on other
necessities such as general welfar, transportations, groceries, etc. We
will examine how rent burden varies by counties in the Charlottesville
and Richmond regions.
ACSData2020 <- st_transform(ACSData2020, 4326)
pal <- colorNumeric("viridis", reverse = TRUE, domain = ACSData2020$MedRentBurdenPercE) # viridis
leaflet(ACSData2020) %>%
#addTiles() %>% # basemap
addProviderTiles("CartoDB.Positron") %>% # CartoDB.Positron, Stamen.Toner, Esri.NatGeoWorldMap
addPolygons(fillColor = ~pal(MedRentBurdenPercE),
fillOpacity = 0.8,
color = "white",
weight = 2,
highlight = highlightOptions(
weight = 3,
fillOpacity = 0.9,
bringToFront = T),
popup = paste0("Region: ", ACSData2020$Region, "<br>",
"County: ", ACSData2020$County, "<br>",
"Tract: ", ACSData2020$NAMELSAD, "<br>",
"Percent: ", round(ACSData2020$MedRentBurdenPercE, 2)))%>%
addLegend("topright", pal = pal, values = ACSData2020$MedRentBurdenPercE,
title = "Percentage of Rent Burden", opacity = 0.7)
ACSData2020 %>%
#filter(Region == "Charlottesville") %>%
ggplot()+
geom_sf(aes(fill=MedRentBurdenPercE))+
scale_fill_viridis_c()+
scale_color_viridis_c()+
ggtitle("Rent Burden in Charlottesville and Richmond, VA", subtitle= "Data from the 2020 American Community Survey (ACS)")+
labs(fill="Rent Burden (%)")
ggsave("Rent Burden in Charlottesville and Richmond.png")
## Saving 7 x 5 in image
The darker the tract the more financially burdened by their rent.
How does the median renters household income relate to the median rent they are paying as proof of rent burden?
ggplot(ACSData2020) +
aes(x = MedRenterHHIncomeE, y = MedRentE) +
geom_point()
## Warning: Removed 32 rows containing missing values (geom_point).
As expected the lower the household income the less the renter can
afford to pay. However does how much they are paying
PovertybyRace <- read_csv("PovertyByRaceData.csv") %>%
mutate(GEOID=as.character(GEOID))
## Rows: 317 Columns: 69
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): NAME, Region, County, FipCodes
## dbl (65): GEOID, TotalPopE, TotalPopM, PovTotalE, PovTotalM, AbovePovTotalE,...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
PovertybyRace
## # A tibble: 317 × 69
## GEOID NAME Region County FipCodes TotalPopE TotalPopM PovTotalE PovTotalM
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 5100301… Cens… Charl… Albem… 003 4639 621 292 177
## 2 5100301… Cens… Charl… Albem… 003 5370 526 218 242
## 3 5100301… Cens… Charl… Albem… 003 3021 400 111 72
## 4 5100301… Cens… Charl… Albem… 003 2011 486 433 226
## 5 5100301… Cens… Charl… Albem… 003 1342 127 12 21
## 6 5100301… Cens… Charl… Albem… 003 5370 677 130 109
## 7 5100301… Cens… Charl… Albem… 003 4255 600 225 108
## 8 5100301… Cens… Charl… Albem… 003 3028 424 104 95
## 9 5100301… Cens… Charl… Albem… 003 2663 438 63 63
## 10 5100301… Cens… Charl… Albem… 003 3189 406 51 47
## # … with 307 more rows, and 60 more variables: AbovePovTotalE <dbl>,
## # AbovePovTotalM <dbl>, IncomeInequalityE <dbl>, IncomeInequalityM <dbl>,
## # WhiteTotalE <dbl>, WhiteTotalM <dbl>, WhiteAboveE <dbl>, WhiteAboveM <dbl>,
## # WhitePovertyTotalE <dbl>, WhitePovertyTotalM <dbl>,
## # WhitePovertyPercE <dbl>, WhitePovertyPercM <dbl>, BlackTotalE <dbl>,
## # BlackTotalM <dbl>, BlackAboveE <dbl>, BlackAboveM <dbl>,
## # BlackPovertyTotalE <dbl>, BlackPovertyTotalM <dbl>, …
TenureByRace <- read_csv("TenureByRaceData.csv") %>%
mutate(GEOID=as.character(GEOID))
## Rows: 317 Columns: 41
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): NAME, Region, County, FipCodes
## dbl (37): GEOID, TotalHHE, TotalHHM, TotalOwnerOcE, TotalOwnerOcM, TotalRent...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
TenureByRace
## # A tibble: 317 × 41
## GEOID NAME Region County FipCodes TotalHHE TotalHHM TotalOwnerOcE
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 51540000201 Census Tr… Charl… Charl… 540 1146 188 270
## 2 51540000202 Census Tr… Charl… Charl… 540 1544 256 190
## 3 51540000302 Census Tr… Charl… Charl… 540 1607 280 633
## 4 51540000401 Census Tr… Charl… Charl… 540 1471 169 721
## 5 51540000402 Census Tr… Charl… Charl… 540 2128 207 1196
## 6 51540000501 Census Tr… Charl… Charl… 540 1666 178 615
## 7 51540000502 Census Tr… Charl… Charl… 540 2033 223 872
## 8 51540000600 Census Tr… Charl… Charl… 540 1275 151 96
## 9 51540000700 Census Tr… Charl… Charl… 540 1651 171 828
## 10 51540000800 Census Tr… Charl… Charl… 540 1451 176 971
## # … with 307 more rows, and 33 more variables: TotalOwnerOcM <dbl>,
## # TotalRenterOcE <dbl>, TotalRenterOcM <dbl>, WhiteOwnersE <dbl>,
## # WhiteOwnersM <dbl>, WhiteRentersE <dbl>, WhiteRentersM <dbl>,
## # PercOfWhtRenters <dbl>, BlackOwnersE <dbl>, BlackOwnersM <dbl>,
## # BlackRentersE <dbl>, BlackRentersM <dbl>, PercOfBlkRenters <dbl>,
## # NativeAmOwnersE <dbl>, NativeAmOwnersM <dbl>, NativeAmRentersE <dbl>,
## # NativeAmRentersM <dbl>, PercOfNaAmRenters <dbl>, AsianOwnersE <dbl>, …
A representation of the percentage of Black renters in and around Charlottesville and Richmond Virginia to compare the regions with the predominately greater than Black population with the predominately rent burdened or exploited, etc.
TenureByRace<-countytracts %>%
left_join(TenureByRace, by="GEOID")
TenureByRace <- st_transform(TenureByRace, 4326)
pal <- colorNumeric("viridis", reverse = FALSE, domain = TenureByRace$PercOfBlkRenters) # viridis
leaflet(TenureByRace) %>%
#addTiles() %>% # basemap
addProviderTiles("CartoDB.Positron") %>% # CartoDB.Positron, Stamen.Toner, Esri.NatGeoWorldMap
addPolygons(fillColor = ~pal(PercOfBlkRenters),
fillOpacity = 0.8,
color = "white",
weight = 2,
highlight = highlightOptions(
weight = 3,
fillOpacity = 0.9,
bringToFront = T),
popup = paste0("Region: ", TenureByRace$Region, "<br>",
"County: ", TenureByRace$County, "<br>",
"Tract: ", TenureByRace$NAMELSAD, "<br>",
"Percentage: ", round(TenureByRace$PercOfBlkRenters, 2)))%>%
addLegend("topright", pal = pal, values = TenureByRace$PercOfBlkRenters,
title = "Percentage of Black Renters", opacity = 0.7)
TenureByRace <- st_transform(TenureByRace, 4326)
pal <- colorNumeric("viridis", reverse = FALSE, domain = TenureByRace$PercOfBlkRenters) # viridis
TenureByRace %>%
filter(PercOfBlkRenters >= 40 ) %>%
#create new dataframe with the above call that datafram in leaflet arguments below
leaflet() %>%
#addTiles() %>% # basemap
addProviderTiles("CartoDB.Positron") %>% # CartoDB.Positron, Stamen.Toner, Esri.NatGeoWorldMap
addPolygons(fillColor = ~pal(PercOfBlkRenters),
fillOpacity = 0.8,
color = "white",
weight = 2,
highlight = highlightOptions(
weight = 3,
fillOpacity = 0.9,
bringToFront = T),
popup = paste0("Region: ", TenureByRace$Region, "<br>",
"County: ", TenureByRace$County, "<br>",
"Tract: ", TenureByRace$NAMELSAD, "<br>",
"Percentage: ", round(TenureByRace$PercOfBlkRenters, 2)))%>%
addLegend("topright", pal = pal, values = TenureByRace$PercOfBlkRenters,
title = "Percentage of Black Renters", opacity = 0.7)
Note the lack of representation of housholds in particular census tracts
#joining the spatial data with the PovertybyRace dataset
PovertybyRace<-countytracts %>%
left_join(PovertybyRace, by="GEOID")
The following is a visual of the Gini Index of Charlottesville and Richmond which entails the index of income inequality in the tracts with the closer a tract is to 1 the closer they are to perfect Inequality and vice versa for 0.
PovertybyRace <- st_transform(PovertybyRace, 4326)
pal <- colorNumeric("viridis", reverse = TRUE, domain = PovertybyRace$IncomeInequalityE) # viridis
leaflet(PovertybyRace) %>%
#addTiles() %>% # basemap
addProviderTiles("CartoDB.Positron") %>% # CartoDB.Positron, Stamen.Toner, Esri.NatGeoWorldMap
addPolygons(fillColor = ~pal(IncomeInequalityE),
fillOpacity = 0.8,
color = "white",
weight = 2,
highlight = highlightOptions(
weight = 3,
fillOpacity = 0.9,
bringToFront = T),
popup = paste0("County: ", PovertybyRace$County, "<br>",
"Tract: ", PovertybyRace$NAMELSAD, "<br>",
"Gini Index: ", round(PovertybyRace$IncomeInequalityE, 2)))%>%
addLegend("topright", pal = pal, values = PovertybyRace$IncomeInequalityE,
title = "Gini Index of Income Inequality", opacity = 0.7)
scatterplot matrix
povertymatrix<- ggpairs(PovertybyRace, mapping=aes(color=Region), columns=c(31, 39, 55, 79), title="Poverty Percentage by Race")
povertymatrix
## Warning: Removed 3 rows containing non-finite values (stat_density).
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 7 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 56 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 21 rows containing missing values
## Warning: Removed 7 rows containing missing values (geom_point).
## Warning: Removed 6 rows containing non-finite values (stat_density).
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 60 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 23 rows containing missing values
## Warning: Removed 56 rows containing missing values (geom_point).
## Warning: Removed 60 rows containing missing values (geom_point).
## Warning: Removed 56 rows containing non-finite values (stat_density).
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 71 rows containing missing values
## Warning: Removed 21 rows containing missing values (geom_point).
## Warning: Removed 23 rows containing missing values (geom_point).
## Warning: Removed 71 rows containing missing values (geom_point).
## Warning: Removed 21 rows containing non-finite values (stat_density).
cvl_rva_measures <- read_csv("cvl_rva_measures.csv")
## Rows: 317 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): county_tract, region
## dbl (2): dissim_wb_block, iso_b_block
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
cvl_rva_measures<-select(cvl_rva_measures,c(1, 4,2,3))
cvl_rva_measures
## # A tibble: 317 × 4
## county_tract region dissim_wb_block iso_b_block
## <chr> <chr> <dbl> <dbl>
## 1 003010100 cvl 0.674 0.145
## 2 003010201 cvl 0.565 0.187
## 3 003010202 cvl 0.485 0.0290
## 4 003010301 cvl 0.210 0.125
## 5 003010302 cvl 0.315 0.0277
## 6 003010303 cvl 0.288 0.0856
## 7 003010401 cvl 0.602 0.278
## 8 003010402 cvl 0.504 0.113
## 9 003010501 cvl 0.356 0.135
## 10 003010502 cvl 0.460 0.281
## # … with 307 more rows
countytracts_updated<-countytracts %>%
mutate(county_tract=paste(COUNTYFP,TRACTCE, sep=""))
cvl_rva_measures<-countytracts_updated %>%
left_join(cvl_rva_measures,countytracts_updated, by="county_tract")
cvl_rva_measures<- cvl_rva_measures %>%
mutate(FipCodes=str_sub(GEOID,3,5)) %>%
mutate(County = case_when(
FipCodes=="540"~"Charlottesville",
FipCodes=="109"~"Louisa",
FipCodes=="125"~"Nelson",
FipCodes=="065"~"Fluvanna",
FipCodes=="003"~"Albemarle",
FipCodes=="079"~"Greene",
FipCodes=="760"~"Richmond City",
FipCodes=="087"~"Henrico County",
FipCodes=="041"~"Chesterfield County",
FipCodes=="670"~"Hopewell City",
FipCodes=="730"~"Petersburg City"
))
cvl_rva_measures<-cvl_rva_measures %>%
mutate(Region=case_when(
FipCodes %in% c("540","109","125", "065", "003", "079")~"Charlottesville",
FipCodes %in% c("760","087","041", "670", "730")~"Richmond"
))
cvl_rva_measures <- st_transform(cvl_rva_measures, 4326)
pal <- colorNumeric("viridis", reverse = FALSE, domain = cvl_rva_measures$dissim_wb_block) # viridis
# cvl_rva_measures %>%
# filter(dissim_wb_block >= 0.5 ) %>%
leaflet(cvl_rva_measures) %>%
#addTiles() %>% # basemap
addProviderTiles("CartoDB.Positron") %>% # CartoDB.Positron, Stamen.Toner, Esri.NatGeoWorldMap
addPolygons(fillColor = ~pal(dissim_wb_block),
fillOpacity = 0.8,
color = "white",
weight = 2,
highlight = highlightOptions(
weight = 3,
fillOpacity = 0.9,
bringToFront = T),
popup = paste0("Region: ", cvl_rva_measures$Region, "<br>",
"County: ", cvl_rva_measures$County, "<br>",
"Tract: ", cvl_rva_measures$NAMELSAD, "<br>",
"Dissimilarity Index ", round(cvl_rva_measures$dissim_wb_block, 2)))%>%
addLegend("topright", pal = pal, values = cvl_rva_measures$dissim_wb_block
,
title = "Dissimilarity Index", opacity = 0.7)
Isolation
#cvl_rva_measures <- st_transform(cvl_rva_measures, 4326)
pal <- colorNumeric("viridis", reverse = FALSE, domain = cvl_rva_measures$iso_b_block
) # viridis
# cvl_rva_measures %>%
# filter(dissim_wb_block >= 0.5 ) %>%
leaflet(cvl_rva_measures) %>%
#addTiles() %>% # basemap
addProviderTiles("CartoDB.Positron") %>% # CartoDB.Positron, Stamen.Toner, Esri.NatGeoWorldMap
addPolygons(fillColor = ~pal(iso_b_block),
fillOpacity = 0.8,
color = "white",
weight = 2,
highlight = highlightOptions(
weight = 3,
fillOpacity = 0.9,
bringToFront = T),
popup = paste0("Region: ", cvl_rva_measures$Region, "<br>",
"County: ", cvl_rva_measures$County, "<br>",
"Tract: ", cvl_rva_measures$NAMELSAD, "<br>",
"Issolation Index ", round(cvl_rva_measures$iso_b_block, 2)))%>%
addLegend("topright", pal = pal, values = cvl_rva_measures$iso_b_block
,
title = "Issolation Index", opacity = 0.7)
#{}
segregationmatrix<- ggpairs(cvl_rva_measures, mapping=aes(color=Region), columns=c(15,16), title="Resedential Segregation")
segregationmatrix
DataforMatrix<-cbind(ACSData2020,TenureByRace, PovertybyRace, cvl_rva_measures)
CvlDataforMatrix<-DataforMatrix %>%
filter(region=="cvl")
RvaDataforMatrix<- DataforMatrix %>%
filter(region=="rva")
rentmatrix<- ggpairs(CvlDataforMatrix, columns=c(21,25,72, 131), mapping=aes(color=region), title="Corelation between Rent Exploitation, Burden, and Resedential Segregation in Black Renters", lower = list(continuous = "smooth"))
rentmatrix
## Warning: Removed 1 rows containing non-finite values (stat_density).
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removing 1 row that contained a missing value
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removing 1 row that contained a missing value
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 2 rows containing missing values
## Warning: Removed 1 rows containing non-finite values (stat_smooth).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing non-finite values (stat_density).
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removing 1 row that contained a missing value
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 2 rows containing missing values
## Warning: Removed 1 rows containing non-finite values (stat_smooth).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing non-finite values (stat_smooth).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removing 1 row that contained a missing value
## Warning: Removed 2 rows containing non-finite values (stat_smooth).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 2 rows containing non-finite values (stat_smooth).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing non-finite values (stat_smooth).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing non-finite values (stat_density).
#upper = list(continuous = "density", combo = "box_no_facet"),
#),, combo = "facetdensity"
#filter(Region ==Charlottesville)
#
# # lower = list(continuous = "smooth"),
# r <- rentmatrix[3,1]
# #
# #
# r +geom_smooth()
# rentma
rentmatrix<- ggpairs(DataforMatrix, columns=c(21,25,72, 67, 115), mapping=aes(color=region, alpha=0.1), title="Correlation between Rent Exploitation, Rent Burden & Select Demographic Groups", lower = list(continuous = "smooth", alpha = 1/40))
#
# r <- rentmatrix[5,1]
# r+scale_color_manual(values=c("#d7191c", "#2c7bb6")
rentmatrix
## Warning: Removed 23 rows containing non-finite values (stat_density).
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 23 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 23 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 23 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 23 rows containing missing values
## Warning: Removed 23 rows containing non-finite values (stat_smooth).
## Warning: Removed 23 rows containing missing values (geom_point).
## Warning: Removed 13 rows containing non-finite values (stat_density).
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 13 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 13 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 13 rows containing missing values
## Warning: Removed 23 rows containing non-finite values (stat_smooth).
## Warning: Removed 23 rows containing missing values (geom_point).
## Warning: Removed 13 rows containing non-finite values (stat_smooth).
## Warning: Removed 13 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing non-finite values (stat_density).
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 4 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 4 rows containing missing values
## Warning: Removed 23 rows containing non-finite values (stat_smooth).
## Warning: Removed 23 rows containing missing values (geom_point).
## Warning: Removed 13 rows containing non-finite values (stat_smooth).
## Warning: Removed 13 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing non-finite values (stat_smooth).
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing non-finite values (stat_density).
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 4 rows containing missing values
## Warning: Removed 23 rows containing non-finite values (stat_smooth).
## Warning: Removed 23 rows containing missing values (geom_point).
## Warning: Removed 13 rows containing non-finite values (stat_smooth).
## Warning: Removed 13 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing non-finite values (stat_smooth).
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing non-finite values (stat_smooth).
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 2 rows containing non-finite values (stat_density).
ggsave("Charlottesville and Richmond Scatter Plot Matrix.png")
## Saving 7 x 5 in image
## Warning: Removed 23 rows containing non-finite values (stat_density).
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 23 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 23 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 23 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 23 rows containing missing values
## Warning: Removed 23 rows containing non-finite values (stat_smooth).
## Warning: Removed 23 rows containing missing values (geom_point).
## Warning: Removed 13 rows containing non-finite values (stat_density).
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 13 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 13 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 13 rows containing missing values
## Warning: Removed 23 rows containing non-finite values (stat_smooth).
## Warning: Removed 23 rows containing missing values (geom_point).
## Warning: Removed 13 rows containing non-finite values (stat_smooth).
## Warning: Removed 13 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing non-finite values (stat_density).
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 4 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 4 rows containing missing values
## Warning: Removed 23 rows containing non-finite values (stat_smooth).
## Warning: Removed 23 rows containing missing values (geom_point).
## Warning: Removed 13 rows containing non-finite values (stat_smooth).
## Warning: Removed 13 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing non-finite values (stat_smooth).
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing non-finite values (stat_density).
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 4 rows containing missing values
## Warning: Removed 23 rows containing non-finite values (stat_smooth).
## Warning: Removed 23 rows containing missing values (geom_point).
## Warning: Removed 13 rows containing non-finite values (stat_smooth).
## Warning: Removed 13 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing non-finite values (stat_smooth).
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing non-finite values (stat_smooth).
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 2 rows containing non-finite values (stat_density).
#ggsave() saves as pdf